################################################################################
#this is a script to create summaries of spatial surfaces and extract data at specific locations

################################################################################
#define the input data directory
solar.dir = "C:/R/course/data/Spatial_data/"

#define the output folder
outfolder = "C:/R/course/data/output/"
################################################################################
#load necessary libraries

#list the libraries needed
necessary=c("adehabitat")
#check if library is installed
installed = necessary %in% installed.packages()
#if library is not installed, install it
if (length(necessary[!installed]) >=1) install.packages(necessary[!installed], dep = T)

#load the libraries
for (lib in necessary) library(lib,character.only=T)

###############################################################################
#EXERCISE
#start working with solar data. In this exercise we want to create a new asc file that
#is the sum of data in a group of asc files

#set the working directory where the asc files are stored
setwd(solar.dir)

#get a list of asc files in the working directory that are to be summed together
solar.list=(dir())

#define some output datasets (i.e. the new summary asc file)
solar_sum = NULL

#Cycle through asc files in the solar.list and create a new asc file "solar_sum"
#that is the sum of the asc files
for (solar in solar.list) {
        #define the path and name of the asc file to be imported
        tfile=paste(solar.dir, solar, sep="")
        #import the asc file
        tasc = import.asc(tfile)
        #now work with the solar data
        #first check if the summary asc file has already been populated with data
        if (length(solar_sum) == 0) { #if there are no values , populate it
          solar_sum = tasc
        } else { #if already has values, add the new values to it
          solar_sum = solar_sum + tasc
        }
        }

#you now have a spatial layer that is the sum of the solar asc files
#this can be seen visually using the image command
image(solar_sum)
#the new "solar_sum" layer can also be exported using the export.asc command
export.asc(x=solar_sum, file=paste(outfolder,"solar_sum.asc",sep=""))

################################################################################
#if you want to extract data at specific locations in the spatial layer...
#first define some locations
long=c(140, 140, 140, 140) #longitude
lat=c(-35, -30, -25, -20) #latitude
pnts=data.frame(long=long, lat=lat) #create a data frame that is combined longitude and latitude
points(pnts, col="blue", pch=16) #plot points on image
#extract the location specific info and append as a new column
pnts$solar_sum = join.asc(pnts,solar_sum)
#write out data into a spreadsheet
write.csv(x=pnts, file=paste(outfolder,"pnts.csv",sep=""))

################################################################################
#if you want to extract data from a grid overlayed on the spatial layer...
#first define the extent of the grid and sampling interval
long=seq(from=140, to=150, by=2) #longitude
lat=seq(from=-40, to=-20, by=2) #latitude
grid.pnts=expand.grid(long=long, lat=lat) #get all possible combinations of latitude and longitude
points(grid.pnts, col="red", pch=16) #plot grid on image
#extract the location specific info and append as a new column
grid.pnts$solar_sum = join.asc(grid.pnts,solar_sum)
#write out data into a spreadsheet
write.csv(x=grid.pnts, file=paste(outfolder,"grid.csv",sep=""))
################################################################################